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Abstract. In this work, a new approach is developed for analysis and-dt?sign of tran- 
sonic airfoils. A set of full-potential-equivalent equations in von Mises coordinates 
is formulated from the Euler equations under the irrotationality and isentropic as- 
sumptions. This set is composed of a main equation for the main variable y, and 
a secondary equation for the secondary variable R. The main equation is solved by 
type-dependent differencing combined with a shock point operator. The secondary 
equation is solved by marching from a non-characteristic boundary. Sample compu- 
tations on NACA 0012 and biconvex airfoils show that, for the analysis problem, the 
present approach achieves good agreement with experimental C p distributions. For 
the design problem, the approach leads to a simple numerical algorithm in which the 
airfoil contour is calculated as a part of the flow field solution. 


7 


: 7 

p-l°- 


1. Introduction 

Transonic flow is a widely encountered phenomenon in aeronautics and astronau- 
tics but is not easy to calculate because the flow field, and the governing equations 
as well, are mixed type. Therefore, transonic computation had little progress until 
1971 when Muraian and Cole developed a type-dependent difference scheme and 
successfully solved the transonic small disturbance (TSD) equation[l]. Since then, 
transonic computation has become one of the most upsurging topics for computa- 
tional fluid dynamicists[2-8]. In 1974, Jameson extended transonic flow computa- 
tion to the full potential (FP) stage by constructing a rotated difference scheme[4]. 
Afterwards, papers were published on transonic computation by solving Euler equa- 
tions^^] and their equivalent streamfunction-vorticity formulation[7,8]. Neverthe- 
less, in spite of the recent active efforts on Euler solvers, the full potential calculation 
is still attractive due to its simplicity, efficiency and sufficient accuracy. 

The von Mises transformation is a type of streamline-based transformation which 
generates a streamwise coordinate system. The von Mises formulation has a number 
of advantages when applied in CFD. For example, one can resolve the problem of 
body-fitting coordinates without performing any grid generation. This is because 
the governing equation (flow physics) and grid generation equation (flow geometry ) 
are combined together in this formulation. Furthermore, the boundary condition on 
the airfoil for the analysis problem is Dirichlet, and a non-iterative design technique 
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can be developed for the inverse problem, leading to simplified numerical algorithms 
and a saving of computer time. Therefore, since Barron[9] connected the von Mises 
transformation with Martin’s approach[10] and solved incompressible 2-D symmet- 
ric flow, numerical simulations based on the von Mises transformation have been 
considerably extended, such as to incompressible liftingfll], axisymmetric[12] and 
design[13] problems, and to transonic flow[14,15]. In addition, Greywall[16] and 
Dulikravich[17] obtained a similar formulation for incompressible and compressible 
flows, respectively. 

However, when extending Barron’s approach[9] to transonic flow, several prob- 
lems appear. For compressible flow, apart from the von Mises variable y, another 
variable, the density />, must be updated in each iteration. But in the transonic 
range, the classical difficulty of double value density-massflux relation still exists. 
Besides, shock waves are not easy to handle in von Mises coordinates either by 
the artificial density technique or by type-dependent differencing. Recently, the 
authors [18] developed a new approach to overcome these difficulties by solving so- 
called full-potential-equivalent equations in von Mises coordinates. The principal 
advances over the previous transonic work[14,15] are as follows: 1) To update den- 
sity, instead of solving the non-linear algebraic Bernoulli equation, a first order 
partial differential equation is solved, thereby avoiding the double density problem; 
2) To handle shock waves properly, a shock point operator in von Mises coordinates 
is proposed and combined with the type-dependent difference scheme so that shock 
waves can be captured correctly; 3) Introducing a concept of generalized density 
linearizes the density equation. 

In the next section, an outline of the mathematical formulation is given. The 
numerical algorithms for analysis and design problems are constructed in sections 3 
and 4. In section 5, sample computations are performed to test the approach, and 
conclusions are given in section 6. 

2 . Flow Equations in Streamwise Coordinates 

Two dimensional, steady, inviscid fluid flows are governed by the Euler equations 
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(2-1) 


where p is density, u and v are velocity components in Cartesian coordinates, p is 
pressure, H = p/p + (u 2 + v 2 )/2 is total enthalpy and 7 is the ratio of specific 
heats, p , u, v and p are norm ah zed by free stream density p^, speed and dynamic 
pressure head Poo^lo w ^l e x and y are scaled by the airfoil chord length. 

Introducing streamfunction Vs such that tp y = pu,i(> z = — pv and substituting 
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T into equation (2 — 1), one gets 


i’l/p + P 

-4’zi’y/P ) + 

4’yH 


-V' z4’y/P 
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= 0. 


( 2 - 2 ) 


Streamfunction ip = V ; (^»l/) can be rewritten in an implicit form, F(x ,y ,4') 0« 

or in an explicit form, y = y{x,ip')- This process is equivalent to introducing von 
Mises transformation: x = <p,y = y{<p,ip)- If the Jacobian J = d(x, y)/d{<f>, ip) ^ 
0, oo, then the transformation is single- valued and (2 - 2) becomes 


1 1 /{py*) + pyri> 
y<t>l(pVv ) 

V H 



(2-3) 


where the total enthalpy H = ^p/p + (1 + y*)/(2p 2 y*). The streamline ordinate 
y, called von Mises variable , is viewed as a function of <f> and ip- The velocity 
components can be easily calculated from u — l/(py^),v = y 4 ,/(pyxp), a f tcr V an< ^ P 
Eire solved. 

It is known that the entropy increase across a shock wave is of third order of 
the shock strength. So, if the shock is not strong, transonic flow can be assumed 
isentropic and irrotational. Replacing the energy equation in (2-3) by the isentropic 
relation and keeping in mind that <p = x , we reduce (2 — 3) to 
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where Moo is free stream Mach number and the last equation is the irrotationality 
condition, = 0, expressed in von Mises coordinates. Substituting (2 - 4c) into 
(2 — 4a) Emd (2 — 46) and expanding (2 — 4d), we get 


— j i) y*y+ ^7 ^ — * 


, 7+1 


pt 


2 p* 


y^yzz - yzVzii> - yzy<p— + Vv> 


2 p Tf 1 p+ 


p ^ M lo P 


= 0, 


(2-5a) 

(2-5b) 


1 


L 


J 



362 

Third International Conference on Inverse Design Concepts and Optimization in Engineering Sciences 
nCIDES-Iin. Editor: G.S. Dulikravich. Washington D.C.. October 23-25. 1991. 

^ ~ 2y z y^,y rt p + (1 + y z )y w — J hy%~~ + J/v(l + y\ ) = 0. (2-5c) 

p p 

Properly manipulating the above three equations can produce several sets of equa- 
tions. Each set has two independent equations for two variables. To make the 
formulation more compact, define generalized density R = p" !+l as an alternative 
to density p. Solving for p z jp and pp/p from (2 - 5a)and (2-56), and substituting 
into (2 — 5c), one gets 


(yp )yxz — 2yrypy t p + (1 + y\)ypp = 0. (2-6a) 

Eliminating y z p from (2 - 5a) and (2 - 56) gives 

y* y\ R t - vA 1 + vl)Rp = (7 + 1 )Ml c ,y zz . (2-cb) 

Equation (2 — 5a) can be rewritten as 


Ml 


yA y p - - y*v% R * = (7 + i)M£,y, 0 . (2-6c) 

Substituting the above y zz ,yz<p into (2 — 5c) and replacing p by R , one obtains 
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It is important to note that (2 66) is linear after introducing the new variable R. 

The term M ^/ R is usually called compressibility factor. 

In principle, any two of the above four equations could be combined as a set of 
equations to solve for y and R. But, in practice, equation (2 — 6a) is always selected 
to solve for y and one of the remaining three equations is selected to solve for R. 
Equation (2 — 6a) is a second order, non-linear, partial differential equation of mixed 
type depending on the local flow property. If the flow is subsonic/supersonic, then 
(2 -6a) is elliptic/hyperbolic. In other words, the mathematical classification of the 
equation is consistent with the physical nature of the local flow. Therefore, (2 - 6a) 
is named the main equation for the corresponding main variable y. Equations (2 - 
66) - (2 - 6d) are called secondary equations for the secondary variable R. Among 
the three secondary equations, (2 66) appears simpler because it is linear and 

hence priority is given to it to accompany the main equation. The main equation 
(2 6a) and one of the secondary equations (2 — 66) — (2 — 6<f) constitute a set 

of so-called full-potential-equivalent equations. They are coupled with each other 
and solved in an alternating and iterative manner. More details and other forms of 
full-potential-equivalent equations can be found in [18]. 
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3. Analysis Problem 

For a symmetric airfoil placed in a transonic stream at zero angle of attack, the 
governing equations (2 - 6a) and (2 — 66) can be rewritten as 

^l!/rz 4* A 2 y T ^ 4" A^y^yp “ 0? (3*1) 

B l R x + B 2 R rl) ^B z (3-2) 

where A 1 = y\, - M^/R,A 2 = -2y z y^,A 3 = 1 + y 2 z ,B 1 = y z y\„B 2 = -y*( 1 + 
y 2 z ),B 3 = ( 74 - 1 )Ml c y zz . The boundary conditions on y are Dirichlet: y = f(x) 


on the airfoil, y = ^ a t infinity, y = 0 on the symmetry line and R = 1 at infinity, 
where f(x) is the airfoil shape function. The computational domain and boundary 
conditions are shown in Fig.l. 

Since the mathematical character of (3 — 1) depends on the local flow property, 
it is necessary to apply Murman and Cole’s type-dependent scheme[l] to solve for 
y. Applying the type-dependent difference scheme to (3 — 1) gives 

A Vij~i + By { j 4- Cyij+x = RHS (3-3) 

where A = p 2 A z - ^(3A 2 ,B = -2/? 2 A 3 + (1 - 3 u)A u C = f3 2 A z + 

RHS = - vAi(yi+ 4- y,_ ij) 4- (1 - i^)A 1 (2y l _i j - y,_ 2 j) 

- ^/3A 2 (t/,+i i; >i - j/i+ij-i - y,*-ij+i + 2 /t — i — i ) / 4 
4* (1 - v)f3A 2 {yi- ij+i - t/i-i j-i)/2, 

and ft = Ax/ A Vq for i — 2,3 , ...,J maz - l,j = 2,3, ..., J mar — 1. The switch 
parameter z/ = 1 for a subsonic point, i/ = 0 for a supersonic point. The resulting 
system of difference equations (3 — 3) has a tridiagonal coefficient matrix so that 
SLOR can be applied by relaxing along vertical lines, sweeping from left to right 
and iterating up to convergence, (see Fig.l) 

After y(x,tfr) is solved from the main equation (3-1) and y*,y0,y r r are properly 
differenced, the secondary equation (3 — 2) can be solved for /2(z,t/>) by marching 
from an initial line other than its characteristic curve. The slope of its characteristic 
curve is difr/dx = —(1 + y^)/(y r y^)- At infinity, d\f>fdz — oo. Thus, left and right 
far field boundaries are characteristic curves and hence cannot serve as initial lines. 
Fortunately, the horizontal boundary is not a characteristic and we can march 
equation (3 — 2) from the top boundary to the airfoil using the condition R = 1 at 
infinity. 

The Crank-Nicolson scheme for (3 — 2) gives 

ARi-ij + BR itj + C/2,+1,; = RHS (3-4) 
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r where RH S = CRi- iJ+1 + BR i<j+1 + AR i+liHl + 4AxB 3 ,A = -B U B = -4 (3B 2 , 
C = Bi,f3 = Ax/ Arp for j = J maz — 1, 3, 2, 1 , i = 2, 3, I maz — 1- The system 
of difference equations (3-4) can be solved row by row from the horizontal far field 
to the airfoil using SLOR, but no iteration is needed because (3 - 2) is linear. After 
R is solved, the pressure coefficient is calculated from 

c - = ihi {RTb - < 3 - 5 > 

However, it has been found after numerical tests that this procedure is efficient 
only for flow in which the shock is weak. For flow with a stronger shock, the itera- 
tions fail to converge. To overcome this difficulty, a special treatment of the shock 
wave is proposed following the ideas of Murman’s shock structure analysis[2,3]. For 
usual transonic flows, the shock wave is approximately normal and the shock jump 
conditions are given by 

1^1 = °> M = 0- (3-6) 

where (...) represents a jump across the shock. Based on this analysis, the difference 
quotient approximations to y tx ,y z 0 at a shock point, i.e. grid point just behind the 
shock, are constructed a s below: 

( yzz)i,j — ^j .2 (ih'+i.i — Vij — otjyi-ij + ftjt/i- 2 ,j) (3- 7a) 


- 4AxA^ 3/,+1,j+1 ~ 1 - 3/..J-1 

- 3y,_i, ;+ i + 3 yi-ij-i + yi-2j+i - yi-2j- 1)- (3-7b) 

where aj is the density jump factor on j th streamline and given by the Rankine- 
Hugoniot relation of a normal shock. (3 - 7a) and (3 - 76) are called shock point 
operator (SPO) in von Mises coordinates. The difference equations (3 — 3) for y 
and (3 — 4) for R are modified using SPO. Numerical experimentation has shown 
that SPO must be applied in the y zz ,yxii> terms of the main equation (3 — 1) and 
in the B 3 term of the secondary equation (3 - 2). SPO is a crucial tool to capture 
shock waves in supercritical transonic flows. 

4. Design Problem 

Similar to the analysis problem, the main equation (2 - 6a) and secondary equa- 
tions (2 66) or (2 6c) can be solved for y and R alternatively: 

My zz + My Z ii> + A 3 y^ = 0, (4-la) 
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B\R Z + B 2 R^, = B 3, (4-lb) 

D\R Z + D 2 R^, =r Ds (4-lc) 


where A x = y% - M% 0 /R,A 2 = -2y z y^,,A 3 = 1 + y 2 z ,B 1 = y z yl„B 2 = -y v (l + 

vl)i B 3 = (7 + l)^yrn-Oi = yAyl~ M L/B),D 2 = -y z y%,D 3 = ( 7 +l)M^y rv ,. 
The boundary conditions are the same as in the analysis problem, except on the 
airfoil, which is unknown. There, the pressure coefficient C pt is specified, hence, 
the generalized density is also specified: 

R t = {\+' 1 Ml 0 C p j2f''+ l)h (4-2) 

On the airfoil surface, the Bernoulli equation in von Mises coordinates leads to 

F(x)yl -y\ = 1 (4-3) 

where 

This is a Neumann boundary condition on the airfoil when solving (4 — la) for y. 
(4 — 2) is a Diricklet boundary condition on the airfiol when solving (4 — lc) for R. 
In addition, on a symmetry line off the airfoil, R^ = 0. 

If streamlines do not intesect each other on the airfoil, then y ^ > 0, and if, fur- 
thermore, F(x) ^ 0 on the airfoil, then equation (4 — 3) gives y ^ = y/(l + y 2 z )/F(x). 
For most practical transonic flows the required conditions are easily satisfied as long 
as C p , is reasonably specified. Differencing y^, we get 

y .,1 = [4y,\2 - y.,3 - 2G(xi)]/3 (4-4) 

where G(x,) = AV>\/[1 + (yz)i,i]/-^'( a: i)- Considering this new boundary condition, 
we modify system (3 — 3) as follows: 

For j = 2, equation (3 — 3) reads Ayii + By{ <2 + Cyi <3 = RHS. Substituting 
(4 — 4) into it, we have 

(■ B + 4i4/3)y, i 2 + (C - A/ 3)y fi3 = RHS + 2i4G(x,)/3 (4-5) 

Replacing the first equation in system (3 — 3) by (4 — 5), solving the resulting system 
and applying (4 — 4), we can obtain the desired airfoil contour /(x,) = y,i without 
further iteration of the airfoil shape. The computational domain and boundary 
conditions are shown in Fig. 2. 

To solve for the secondary variable R, two secondary equations (4 — 16) and 
(4 — lc) are available. For equation (4 — 16), the marching procedure is the same 
as in the analysis problem, while for equation (4 — lc), the marching procedure is 
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I different. The slope of its characteristic curve is dip/dx = — {y z y^)l{y\, — M^./R). 
At infinity, d^’/dx = 0. So, the horizontal far field boundary is a characteristic 
curve, but the vertical boundaries are not. Therefore, the marching process can be 
carried out from left to right. 

Crank- Nicolson scheme for (4 — lc) gives 

ARij-i + BRi,j + CRij+i = RH S ( 4 - 6 ) 

where RH S = CRi-ij-y + BRi-ij + AR,_ 1>J+1 + 4A xD 3 ,A = -f3D 2 ,B = 4D l5 
C = (3D 2 ,(3 = Ax/Aif), for i = 2,3 J ma * - 1 ,j = 2,3,..., J max - 1. 

For the first equation in system (4 - 6), the boundary conditions R tl = R,(x t ) 
on the airfoil and i?, i = R{ 2 on symmetry line should be imposed. It is noted that 
y z in D 2 should be type-dependent differenced with SPO to keep consistency with 
the main equation. 

Both (4 — 16) and (4 — lc) have been coupled with (4 — la). Numerical experiments 
have shown that (4 — lc) gives better accuracy than (4 — 16). This is reasonable 
because the boundary condition on the airfoil is considered not only in the main 
equation (4 — la), but also in the secondary equation (4 — lc), while it is not suitably 
considered in the secondary equation (4 — 16). However, the price to pay is more 
iterations because (4 — lc) is non-linear. 

5. Sample Computations 

The approach developed here is applied to calculated transonic flows for both 
analysis and design problems. Only symmetric airfoils at zero angle of attack are 
considered, but both subcritical and supercritical Mach numbers are included. In 
the computational domain, a 65x33 uniform mesh covers — 2 < x < 3, 0 < if> < 2.5 
and the airfoil is placed between 0 and 1. For higher Mach numbers, a 80x31 mesh 
has been used. The computational domain and boundary conditions are shown in 
Figures 1 and 2. 

Figures 3 and 4 are comparisons of calculated C p distributions of NACA 0012 with 
experimental data at NAE[19] for M 0 0 = 0.490 and at ONERA[19] for M = 0.803. 
Figure 5 indicates the calculated C p distribution of a 6 percent biconvex airfoil at 
Moo = 0.909 compared with experimental data at NASA[20). From these plots we 
can see that the present approach is able to accurately predict C p distributions on 
airfoils in transonic flows. The agreement between computed pressure and available 
experimental data is quite satisfactory. For supercritical transonic flows, the shock 
wave can be captured by the presently proposed type-dependent scheme with SPO. 

Figure 6 shows the designed contour of a 6 percent biconvex airfoil compared 
with the exact shape[21]. The specified C p distribution on the airfoil comes from 
experiments at NASA[20] for = 0.909. Figures 7 and 8 give designed NACA 
0012 contours compared with the exact shape[21]. The specified C p is from NAE[19] 
for Mqo = 0.490 and ONERA[19] for = 0.803. Here, we can see that the present 
approach is capable of designing airfoil contours with satisfactory accuracy. 
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6. Conclusions 

1) The newly developed approach based on the full-potential-cqui valent equations 
in von Mises coordinates is able to solve transonic flows for both analysis and design 
problems. 

2) The fuil-potential-equivalent equations are composed of a main equation for 
the corresponding main variable, streamline ordinate y, and a secondary equation 
for the related secondary variable, generalized density R. 

3) The type-dependent difference scheme with shock point operator is effective 
to solve the main equation for y and the shock point operator is crucial to capture 
shock waves in supercritical transonic flows. 

4) The secondary equation can be solved for R by marching from a certain non- 
characteristic, density-specified boundary. Crank-Nicolson scheme proves to be use- 
ful to march such a equation. 

5) For analysis problems, the boundary condition on the airfoil is Dirichlet, which 
is easy to implement. 

6) For design problems, the airfoil contour can be obtained in a non-iterative 
manner because it is a part of the solution of the main equation. 


References 

1. E. M. Murman and J. D. Cole, Calculation of Plane Steady Transonic Flows, 
AIAA J., Vol.9, 114-121 (1971) 

2. E. M. Murman and J. A. Krupp, Solution of the Transonic Potential Equation 
Using a Mixed Finite Difference System, Lecture Notes in Physics, Vol.8, ppl99-206, 
Springer- Verlag, Berlin, 1974 

3. E. M. Murman, Analysis of Embedded Shock Waves Calculated by Relaxation 
Methods, AIAA J., Vol.12, 626-633 (1974) 

4. A. Jameson, Iterative Solution of Transonic Flows over Airfoils and Wings, 
Including Flows at Mach 1, Communications on Pure and Applied Mathematics, 
Vol.27, 283-309 (1974) 

5. R. M. Beam and R. F. Warming, An Implicit Finite Differnce Algorithm for 
Hyperbolic System in Conservation Law Form, J.of Comp. Phys., Vol.22, 87-110 
(1976) 

6. J. C. Steger, Implicit Finite Difference Simulation of Flow around Arbitrary 
Two Dimensional Geometries, AIAA J., Vol.16, 676-686 (1978) 

7. M. Hafez and D. Lovell, Numerical Solution of Transonic Stream Function 
Equation, AIAA J., Vol.21, 327-335 (1983) 

8. H. L. Atkins and H. A. Hassan, A New Stream Function Formulation for the 
Steady Euler Equations, AIAA J., Vol.23, 701-706 (1985) 

9. R. M. Barron, Computation of Incompressible Potential Flow Using von Mises 
Coordinates, J. of Math, and Comp, in Simulation, Vol.31, 177-188 (1989) 


L 


J 



368 


Third International Conference on Inverse Design Concepts and Optimization in Engineering Sciences 

flCIDES-Iin. Editor: , C .S. Puli ta avi ch . W ashington P.C» October 23-21 199L 

r i 

10. M. H. Martin, The Flow of a Viscous Fluid I, Archives for Rational Mechanics 
and Analysis, Vol.41, 266-286 (1971) 

11. R. K. Naeem and R. M. Barron, Lifting Airfoil Calculations Using von Mises 
Variables, Communications in Applied Numerical Methods, Vol. 5, 203-210(1989) 

12. R. M. Barron, S. Zhang, A. Chandna and N. Rudraiah, Axisymmetric Po- 
tential Flow Calculations, Part 1: Analysis Mode, Communications in Applied 
Numerical Methods, Vol. 6, 437-445(1990) 

13. R. M. Barron, A Non-Iterative Technique for Design of Aerofoils in Incom- 
pressible Potential Flow, Communications in Applied Numerical Methods Vol 6 
557-564(1990) 

14. R. M. Barron and R. K. Naeem, Numerical Solution of Transonic Flows on 
a Streamfunction Coordinate System, Inti. J. for Num. Methods in Fluids Vol 9 
1183-1193 (1989) 

15. R. K. Naeem and R. M. Barron, Transonic Computations on a Natural Grid 
AIAA J., Vol.28, 1836-1838(1990) 

16. M. S. Grey wall, Streamwise Computation of 2-D Incompressible Potential 
Flows, J. of Comp. Phys., Vol. 59, 224-231 (1985) 

17. G. S. Dulikravich, A Stream-Function-Coordinate (SFC) Concept in Aero- 
dynamic Shape Design, AGARD VKI Lecture Series, May 14-18, 1990 

18. C.-F. An and R. M. Barron, Numerical Solution of Transonic Full-Potential- 
Equivalent Equations in von Mises Coordinates, to appear 

19. J. J. Thibert, M. Grandjacques and L. H. Ohman, Experimental Data Base 
for Computer Program Assessment, AGARD AR-138, pp.Al-1 — Al-36, 1979 

20. E. D. Knetchtel, Experimental Investigation at Transonic Speeds of Pres- 
sure Distributions over Wedge and Circular-Arc- Airfoil Sections and Evaluation of 
Perforated Wall Interferance’, NASA TN D-15, 1959 

21. I. Abbott and A. E. von Doenhoff, Theory of Airfoil Section, Dover Publ. 
Inc., N.Y., 1959 


L 


J 











